library(dplyr)
library(ggplot2)
library(rriskDistributions)
library(haven)
generate_non_normal_beta_dist = function(n,q_vect,min,max){
  # y = a + bx, when x=min,y=0; when x=max,y=1
  # x = (y-a)/b
  b = 1/(max-min)
  a = -b*min
  q_vect_normalized = q_vect * b + a
  rst = get.beta.par(q=q_vect_normalized)
  beta_dist = rbeta(n,rst[1],rst[2])
  beta_dist_non_normalized = (beta_dist-a)/b
  return(beta_dist_non_normalized)
}

v1 = generate_non_normal_beta_dist(100000,c(2,17,18),0,20)
v2 = rnorm(100000,mean=mean(v1),sd=sd(v1))

sample_twice_from_same_pop = function(pop,sample_size_1,sample_size_2){
  v1 = sample(pop,sample_size_1)
  v2 = sample(pop,sample_size_2)
  return(list(v1=v1,v2=v2))
}

sample_and_t_test_trial = function(pop,sample_size_1,sample_size_2){
  rst=sample_twice_from_same_pop(pop,sample_size_1,sample_size_2)
  t.test_rst=t.test(rst[['v1']],rst[['v2']])
  return(t.test_rst$p.value)
}
p_rst_vect = c()
p_rst_vect2 = c()
for(i in 1:100000){
   p_rst_vect[i] =sample_and_t_test_trial(v1,20,25)
   p_rst_vect2[i] = sample_and_t_test_trial(v2,20,25)
}